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ABSTRACT 


In-Situ measurements of the speed of sound in the upper ocean have 
revealed the existence of significant dispersion and large fluctuations 
over the frequency range 25-80 kHz. The near-surface values of c ranged 
from +6 M/sec to -3 M/sec relative to the bubble-free value, with a 
maximum estimated error of 0.5 M/sec. It was possible to identify 
bubbles of "surface" radius centered around 54 microns down to 4.3 meter 
depth as well as a population centered around 124 microns (at 4.3 M) 
found at all depths. The speed fluctuation showed near-Gaussian 
probability density functions except at the dispersion center frequen¬ 
cies. The standard deviation of the speed varied from 0.27 M/sec for 
58.0 kHz to 0.52 M/sec at 69.6 kHz. The phase remained temporally 
correlated only for 1.46 sec., whfch was very close to the correlation 
times of the oceanographic variables. The power spectral density of 
the varying phase showed its strongest values at frequencies less than 
0,5 Hz.; this was presumably due to bubbles entrained at orbital frequen¬ 
cies associated with the surface wave system which had maximum energy 
in the same frequency range. 
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I. xINTRODUCTION 


The objective of this research was to take in-situ acoustical 
measurements in the surface layer of the ocean in order to determine 
the existence of dispersion and fluctuations of the speed of sound. 

The anticipated source of frequency-dependent speed of sound in sea-water 
are bubbles, of various radii, of numbers dependent on the depth. A 
knowledge of this dispersive behaviour is of importance for sound 
propagation in near-surface ducts where a change in the speed of sound 
with depth would result in.a change in the refraction of the sound ray. 

If bubbles exist near the sea surface a sound ray would bounce more 
often from the surface over a constant horizontal distance as compared 
with the bubble-free case. 

The problem has already been treated in theory as in the analysis 
worked out by Meyer and Skudrzyk (1950). Fox, Curley and Larson have 
examined in the laboratory environment the phase velocity in water 
containing a high density of air bubbles (1954). Glotov et al. (1962) 
studied a population of bubbles generated by wind-agitated breaking 
waves in a laboratory tank. Of the few known data actually taken in 
the ocean those given by Buxcey et al. (1965) and analyzed in the open 
literature by Medwin (1970) will be referenced here. 

Since virtually no in-situ ocean measurements of bubble effects were 
available, this research was proposed and carried out using the U. S. 

Navy Undersea Research and Development Center (NUC) Oceanographic 
Research Tower in Mission Bay/San Diego as a test platform. 

A stable, accurately known, continuous wave acoustical signal was 
received by two spatially separated hydrophones mounted on the sound 





axis and their output was compared in phase by means of a digital phase 
meter. The phase difference between the two hydrophones was recorded 
at each frequency. The frequency was incremented in steps of approxi¬ 
mately two or four kHz over a range from 25 to 80 kHz and for four 
different water depths. The average results have been interpreted in 
terms of a change in the speed of sound with frequency at a constant 
depth, and with depth at a constant frequency. The statistics of the 
phase change with time have also been calculated. 

At the same time, and at the same location, the local sound amplitude 
modulation due to the spatial and temporal changes in the medium has 
been examined (thesis by LCDR W. J. Smith, December 1971). In addition, 
field measurements were taken to determine the fluctuating temperature, 
salinity, turbulent velocity in x, y, and z-directions, speed of sound 
as given by a velocimeter and wave height in order to allow a better 
judgment of the acoustic-oceanic interaction (theses by Lt. M. Bordy, 

LCDR C. Duchock and Lt. H. Seymour, March 1972), 



II. INSTRUMENTATION 


A steel-pipe frame shaped like a "square C" of height and length 
six ft was used to mount the transducer and the hydrophones as shown 
in Fig. 1. To prevent the frame from acting like a huge tuning fork 
and thereby displacing the hydrophones, the frame was put under tension 
halfway down and at the open end by spring-loaded wires. To minimize 
the reflections from the frame caused by side-lobes of the transducer 
the upper and lower members were covered with acoustic absorbent rubber 
(SoAB) at the specular points. 

Opposite to the open end of the frame and at the center an USRD F27 
unidirectional transducer was mounted. At the sound axis of this trans¬ 
ducer an Atlantic Research LC-10 hydrophone was placed 78 cm from the 
first one and a total of 173 cm from the F27. This distance between the 
hydrophones v/as chosen to make the measuring path long enough to measure 
the effect of the bubbles while at the same time limiting the overall 
size of the frame. 

The transmitting part of the instrumentation (Fig. 2) consisted of 
the highly stable General Radio Coherent Decade Frequency Synthesizer 
Type 1162-A generating a sinusoidal voltage of two volts rms and variabl 
in frequency. This signal was then amplified by a Hewlett Packard Power 
467A Amplifier to 14,1 V rms and then impressed across the transducer 
emitting the sound wave. The receiving part of the instrumentation 
(see Fig. 2) consisted of hydrophone #1 with a sensitivity of -125 dB 
re 1 volt/microbar and hydrophone #2 with a sensitivity of -106 dB re 
1 volt/microbar at 40 kHz which picked up the acoustical signal. The 
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Fig. 2. Transmitting and Receiving Instrumentation 































two signals were amplified by 30 dB by an NUS Pre-Amplifier Model 
2010-030 before being sent through (150) ft of waterproof shielded 
cable. Both signals were then bandpass-filtered by a digitally tuned 
variable Krohn-Hite Model 3322 Filter. The received signal was always 
bandpass-filtered with the center frequency equal to the frequency of 
the transmitted signal through a passband of i 300 Hz, with -48 dB per 
octave attenuation on either side of the passband. After the Bandpass- 
Filter the signals were fed into a Wiltron Digital Phase Meter Model 355 
which measures the phase angle between the two ac voltages of the same 
frequency and provides a visual display as well as a dc voltage output 
proportional to the measured phase angle at 10 mV per degree. This dc 
voltage was then FM recorded on a Sangamo Model 3500 fourteen track 
magnetic tape recorder simultaneously with the oceanographic data at 
a recording speed of 1 7/8 inch pe.r sec. The dynamic range of the 
recorder utilized was .01 to 2,5 volts, with all the readings falling 
within this range. 
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III. OCEANOGRAPHIC COMPONENTS 


It was desired to define the volume of water carrying the acoustical 
signal in terms of the measurable parameters temperature, salinity, 
particle velocity, speed of sound as given by a velocimeter, and surface 
wave height. Therefore ocean sensors were mounted close to the acoustic 
system but not interfering with its performance (see Fig. 3 for geometry). 
Three thermistors (T1, T2, T3) measured the temporal variation in the 
temperature, and the temperature probe CT4) of the Bissett-Berman 
apparatus gave absolute readings between 14° and 19° C. The Bissett- 
Berman salinometer yielded the temporal variation in the salinity (S) 
for values ranging from 37.5 to 39.5 ppt. The speed of sound c = 
c( T, S, depth) was taken by a Ramsey SVTD probe that provided a temporal 
variHion around a DC-value. The particale velocity was obtained with 
an electromagnetic flowmeter reading one horizontal component (u) or (v) 
and the vertical component (w). A pressure wave height indicator and 
a Baylor gauge gave information on the instantaneous wave height. 
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IV. EXPERIMENTAL SETUP 


A continuous wave method was used instead of a pulse-echo system 
in order to avoid possible errors related to pulse rise time, definition 
of the leading edge of the pulse, distortion of the pulse when passing 
through the medium because of frequency-dependent attenuation effects 
and the difficulty to determine the phase shift on reflection (see Ref. 
7). This concern was particularly true since we wanted measurements 
over a wide frequency range. 

The method used to measure the speed of sound is actually derived 
from a phase interference approach. The frequency of the signal applied 
to the transducer is adjusted by means of the frequency synthesizer 
until a wanted phase angle is approached at the output of the digital 
phase meter. Since the Sangamo recorder required a signal of at least 
100 mV and at the same time the digital phase meter is most accurate in 
its reading for phase angles close to zero, the frequency was always 
adjusted to get a phase value of approximately 20 degrees. Since the 
digital phase meter compare?, the signal at input B with respect to the 
signal at input A and cannot detect integer number of wave lengths 
difference between the two signals, the displayed phase determines the 
fractional wave length by which signal B leads or lags signal A depending 
on the sign. This means that an integer number of wave lengths plus 
this fractional part just fits into the separation distance between the 
two hydrophones. Further increasing the frequency will finally change 
the phase by 360 degrees or step up the integer number of wave lengths 
by one and thus produce the next frequency to be recorded. 




Because bubble populations were expected to be large in the radius 
range 40-150 microns we examined the frequency range from 25-80 kHz. 

The 78. cm separation between the two hydrophones caused the number 
of wave lengths within this spacing to be incremented by one approxi¬ 
mately every two kHz. In order to exclude possible errors in the results 
due to temporal changes in the medium it was necessary to examine the 
total frequency range in as short a period of time as possible. There¬ 
fore it was decided to accept 20 minutes as an upper limit for 
stationarity of the medium. The time was also the lower limit to find 
the frequency spectrum of the ocean wave height and turbulent velocity. 

This time allowed us to examine the frequency range in steps of four 
kHz. and to record the time-variable phase difference between the hydro¬ 
phones for one minute at those frequencies. 

In addition to examining the speed of sound with changing signal 
frequency at a constant depth, readings were taken of the variation of 
the speed for sound with varying depths at a constant frequency. Four 
different depths were chosen in order to observe near-surface to near¬ 
bottom effects: 4,3, 9.3, 14.3, 7.3 m. These depths were chosen to 
minimize interference from the underwater beam structure of the tower. 

The order was picked to decrease and to check on the possibility of 
longtime variations in the medium. The depths are in meters below mean 
water, where mean water was 6.7 m below the concrete deck of the NUC tower. 

All the readings were taken within 4 1/2 consecutive hours at night 
and early in the morning on October 22, 1971. The obtained data are 
given in Appendix A. 
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V. DATA REDUCTION 


Five twenty-minute runs (split into ten two-minute individual 
records) make up the recorded analog data. The approach used was to 
read this information out on a strip-chart recorder, to digitize the 
printed analog data, and to use a computer to carry out a statistical 
and spectral analysis. 

The phase fluctuation and the other oceanographic parameters were 
initially recorded at 1 7/8 inches per second. The tape recorded data 
were played back at 7 1/2 inches per second and first recorded on an 
eight channel Brush Mark 200 strip-chart recorder in order to allow 
visual comparison between the various sensor outputs (Fig. 4). 

The tape recorded data was then played back a second time at 7 1/2 
inches and recorded on a two channel Brush strip-chart recorder. This 
compressed the 20 minute run to five minutes and resulted in an increase 
in frequency by a factor of four. The strip-chart recorder was set at 
25 mm/sec. chart speed and 10 mV/div sensitivity, corresponding to 
1 degree/di V. 

The strip-chart records for each frequency at each depth were digitized 
using the Fleet Numerical Weather Central facility tracing digitizer at 
Point Pinos, California. The digitizer sampled the trace in increments 
of hundredths of an inch in the y-direction for each one one-hundredth 
of an inch advance in the x-direction and recorded this on a nine-track 
magnetic tape. Every record was digitized by first following the zero 
Volt reference line on the strip-chart and then moving the tracer 
vertica-Hy at the start of a record from the base-line onto the trace. 
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By doing this it was assured that the digitized data points were 
absolute values referenced to zero Volt. Thus, the scaling factor in 
the y-direction is 

.254 mm/data point = .254 degree/data point 

resulting in 

.254 mm/d,p. X ^ sec/mm = .01016 brush sec/d.p. 

= ,04064 real sec./data point 

This implies that the sampling rate is 25 Hz and the Nyquist rate is 
then 12.5 Hz. 

The digitized data were read out from the magnetic tape and punched 
on data cards using the Fleet Numerical Weather Central's CDC 5000 
computer (see Computer Program 1). The data cards could then be used 
for analysis on the USNPGS IBM 360 digital computer. 



VI. COMPUTATIONAL METHOD 


Since we were working with a continuous wave method the signal 
received by the first hydrophone can be written as 

y^ = a sin (ojt) 

and the signal received by the second hydrophone B has then to be written 
yg = b sin (cot-kx) 

o .p 

where k = the wave number = ^ 

c c 

f = the frequency of the applied signal 
c = the phase wave velocity 
X = the separation between the hydrophones 

Since we compare the two signals with respect to phase by using the 
digital phase meter we have to find values for co for which 
kx = ij; 

where if/ is the total phase difference between the two signals due to 
an integral number of wave lengths and a fraction of the wave length. 

We write 

kx = C2iT)n + (j) n = an integer 

kx = X = (2TT)n + (f) 

tox 

(.2 tt )n + ()) 


c 


fx 



where <|) is in radians 









The final result will be 


c 


fx 


n + 




360 


where (p is in degrees 


Therefore the speed of sound may be obtained by using the measured 
frequency at which a phase difference <|) was noted. In addition, x and 
n must be known. The integer n can be found by calculating c from the 
Wilson equation CKinsler and Frey, Eq. 15.1) using readings for the 
temperature and the salinity. We find 


and then 
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VII. BUBBLE THEORY 


The speed of sound through a fluid is given by 


where w = sound frequency 

k = Btt/x = propagation constant in bubble-free water 

Following the analysis carried out by Morse and Feshbach, an incident 
acoustic pressure, of wave length very large compared to the bubble 
radius R, having ^ velocity potential , will produce a total 

field averaged over all possible sizes and configurations of bubbles 
when written in differential equation formulation of 

< 4* > + k^^^(r) <4^> = 0 


where 

k^Cr) 


k + 



n(r,R) R dR 

(Wq/w)^ - 1 - i5(a)^/w) 


r = coordinate position of scattering region 
k|^ = propagation constant in bubbly water 
R = bubble radius 

= resonant frequency of bubble 
$ = damping constant of bubble 

n(R)dR = number of bubbles per unit volume in radius increment 
R to R + dR 

Ttie quantity — ^ ' - is the index of refraction in the bubbly region 
and the real part defines the local speed of sound. This term is 


Ztt 

k 



R n(R) - 1] dR 

- 1]^ + 


Reikj^Cr)} 


k + 
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c 


k 


k + 



R n(R) [( 0 )^ 0 ) - 1] dR 

I(Wq/w)^ - 1]^ +6^(a)^/u))^ 


1 



£ R^- nCR^OE^Q/f^i)^ - 1] dR. 

i [(“q/^^)^ “ 1]^ + 


where 

nCR) dR 


nCR^)dR. 


ti(R) dR 

4/3 ttR^ 

a(R^.) dR. 
4/3 ttR.^ 


uCR) dR = fractional air to water volume for the increment between 
R and R + dR 

uCR.) dR. = fractional air to water volume for the increment between 
^ ' R^. and R^ + 1/n 


Therefore we can write 


c 




1 


uCR^.) dR. 



-1] _ 

EC^o./'*^)^ - 1]^ + 


This is the dispersion relation that predicts a decrease of sound speed 
caused by bubbles of resonance frequency greater than the incident 
frequency and an increase of sound speed for bubbles with lower resonant 
frequencies than the impinging sound. 













The resonance frequency of clean air bubbles In water is given by 
(Eckart, 1945) 

^0 " ~ (V2TTRQ)I(3YpQ/pQ)(g/a)J'^^'''^ 

where y = ratio of specific heats of bubble gas (=1.402 for air) 

= ambient pressure at bubble depth 
= density of surrounding water 

The factor g is due to the influence of surface tension, it differs from 
unity only for very small bubbles, for the smallest bubbles of this 
study (R^ = 50 microns) g = 1.02. The factor modifies the value 
of Y, i.e., is the effective ratio of specific heats and lies for 
the smallest bubbles of this study in the range 1,24 < ya"^ ^1.40. 

The two effects are to some extent mutually counteracting in their 
effect on the resonance freqcjercy and therefore the term yg/a will 
be assumed to equal 1,33. 

The damping factor as a function of bubble resonance frequency 

as given by Devin will be approximated by 

0,322 

6. = ,027 X 10^^ 
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VIII. EXPERIMENTAL RESULTS 


Five different freqency sweeps have been carried out, runs 6 and 7 
at a depth of 4,3 m and the remaining three, runs 8, 9, and 10, at 
9.3, 14.3 and 7.3 m depth, respectively. 

Run 6 covered a total frequency range from 25 - 80 kHz in steps of 
approximately two kHz. to be able to detect the detailed variations 
in the speed of sound; for the remaining runs readings were taken every 
four kHz, in order to examine the frequencies of interest in as short 
a period of time as possible. Run 7 started on hour after the end of 
run 5, but for the five runs reported here all the readings were taken 
within 4 1/2 consecutive hours between 0350 and 0820 on October 22, 1971. 

Single temperature and salinity readings taken at each depth and 
a single bathythermograph record taken at 8 a.m. (Fig. 5) allow one 
to calculate the speed of sound for bubble-free water; the Wilson 
equation is used in the following form: 


c = 1449 + 4.6t - 0.055t^ +0.0003t^ + (1 .39-,0l2t)(S-35) + 0.017d 


where c = speed of sound in meters per second 
t = temperature in degrees centigrade 
S = salinity in parts per thousand 
d = depth in meters 

The bubble-free values obtained for 4.3, 9,3 and 14.3 m are used as 
references for the dispersion curves and are included in Fig. 6, 

The calculated values for the speed of sound based on the obtained 
data are examined and analyzed with respect to a frequency — and/or depth 
dependence. If there is any change in c with frequency then this change 
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is attributed in general to either molecular relaxation or bubbles 
since for the sea«water medium no other frequencyndependent mechanisms 
are known at the frequencies of this study. But the magnitude of the 
molecular relaxation effect can be shown to be negligible (a frequency 
shift from 20 to 70 kHz. results in a Ac2<0.3 m/sec.) as compared to 
the bubble effects. 

A, FREQUENCY BEHAVIOUR 

All the calculated values of c as a function of frequency for the 
different runs have been plotted on Fig. 6 for comparison. They should 
be interpreted by making use of the reference value for the non-bubbly 
medium which, because of the temperature change, shifts by as much as 
4.2 m/sec. in going from 4.3 m to a depth of 14.3 m. 

The first point of interest is the dispersive variation of c at a 
single depth. The first run (6) was carried out at a depth of 4.3 m 
and examined frequencies from 30.7 up to 79.2 kHz during a 40 minute 
period of time. Most significant is the strong increase in the speed 
of sound by as much as 6 m/sec. compared to the non-bubbly value for 
c of 1510.47 m/sec, and caused by the peak at 34.7 kHz followed by a 
smaller negative change of 3 m/sec at 40.6 kHz. For frequencies from 
47 to 60 kHz the speed of sound shows nearly no dispersive behaviour. 
The higher frequencies are characterized by another peak in the c-curve 
roughly one fourth the strength of the peak at 34.7 kHz. 

Run seven is interesting since it covered the same frequency range 
and the same depth as run six but for the major part now in steps of 
four kHz rather than 2 kHz. Since this allows a comparison between the 
two runs, a reminder is necessary that there is a time lag of one hour 
between the end of run six and the start of run seven. This time delay 


) 




may account for the slight shift in the dispersion curve downward and 
to the right. The shape matches the previous one pretty well, with two 
peaks clearly discernible. 

In run eight the depth is increased by five meters as compared with 
the previous run. The main feature is a resonance pattern with a weaker 
maximum at 36.8 kHz rather than 34.7 kHz at shallower depth. 

Run nine, at the greatest depth (14.3 m), differs in the pattern 
for the speed of sound over the frequency range by showing an increase 
in c over a wide frequency band with a maximum at approximately 41 kHz 
and no second peak. The obvious shift in the level is caused by a 
decrease of the non-dispersive value of c from 1510.47 at 4.3 m to 
1506.27 m/sec at 14.3 m. 

Run ten was made back at a depth of 7.3 m. The two readings at 
57,7 and 61.59 kHz are believed to be non-representative since a boat 
approaching the tower injected a large number of bubbles into the volume 
under consideration. Excluding these data points, the dispersion curve 
shows a main peak around 35.5 kHz. 

Summarizing the dispersion data it can be stated that: 

1) The change of the speed of sound with frequency at sea was observed 
to be significantly stronger than expected. 

2) The change in the speed of sound does not follow the usual dispersive 
pattern for a single dominant bubble regime. 

3) The speed of sound curve at 4,3 m has two peaks. The frequency of 
the lower peak changes from 34.7 kHz at 4,3 m to 41 kHz at 14.3 m, which 
is a shift toward higher frequencies as expected for surface-generated 
bubbles carried downwards "isothermally" by turbulence. However, the 
shift is smaller than expected: isothermally convected identifiable. 




non-diffusing, surface-initiated bubbles which resonate at 35 kHz at 
4.3 m would resonate at 54 kHz and not at the observed 41 kHz when 
they reach 14.3 m. There are two explanations that can be proposed 
for this discrepancy: 

a) the presence of the higher frequency peak in the near-surface 
measurement has shifted the low frequency peak to higher frequencies 
than it would have had by itself. The absence of the second peak in 
the deeper measurement (due to more rapid disappearance of smaller 
bubbles) permits the low frequency peak to more accurately identify 
the resonance frequency of these bubbles. We herein assume that near 
the surface v/e encountered a bubble population made up of larger surface¬ 
generated bubbles and smaller bubbles introduced primarily by continental 
aerosols that drop into the water. To check on the dispersion curve that 
results from a specific bubble distribution, a bubble population was 
created peaking around 60 and 120 microns as given in Fig, 7. Subse¬ 
quently the computer program (3) was written calculating the speed of 
sound dispersion based on the equations derived in Section VII and using 
the above mentioned volume concentration of bubbles. The result is a 
dispersion curve given in Fig. 8 which is in good agreement with our 
data obtained at 4.3 m. The observed dispersion curves for the other 
depths support this hypothesis, since only part of the larger bubbles 
are carried down, decreasing in number with increasing depth whereas the 
smaller bubbles rapidly disappear due to gas diffusion out of the 

bubbles. 

b) the larger bubbles originate in the volume or at the bottom, rather 
than at the surface. As they rise they lose gas by diffusion, so that 
when they reach the surface they are not as large as predicted by the 
non-diffusion assumption and are thereby at higher frequencies than expected. 
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Speed of Sound Dispersion for Bubble Population of Fig. 




















B. DEPTH BEHAVIOUR 

Next the variation in c with depth for a constant frequency will 
be considered. On the assumption that the bubble densities were not 
a function of time during the 2 1/2 hours of runs 7~10, all these data 
have been rearranged to display the variation in the speed of sound 
with depth at various constant frequencies as given in Table I. These 
values have been used to plot the speed of sound versus depth in Fig, 9. 
These straight line segments reflect the major variations in c with 
frequency as described in the last section. 

In order to examine the variation of the speed of sound versus 
depth, it is necessary to subtract out first all possible temperature, 
salinity, and pressure-based changes with depth. To do this, use is 
made of the following coefficients for variation with 

temperature aT ~ ^ ^ ft/(sec) C°F) 

salinity fs" " ^ ^ ft/(sec) (ppt) 

depth aI " ft/Csec) (ft) 

The salinity changed by -.1 ppt over 10 m and therefore accounts for a 
change of c with S of -.4 meter per second. It is necessary to add to 
this a change of +.17 m/sec due to the depth effect such that the total 
effect of salinity and depth on c will be -.23 m/sec when going from 
4,3 to 14.3 m; this small effect will be neglected. 

The temperature-caused changes in the speed of sound have to be 
taken into consideration. The BT curve sketched in Fig. 5 is approxi¬ 
mated by straight line elements and is used with the Wilson equation 
to calculate values for the speed of sound at different depths for the 
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TABLE I. Speed of Sound in m/sec 


-^,Depth^^(^ 
Freq. (kHz) 

4.3 

4,3 

7.3 

9.3 

14.3 

28.97 




1505.73 

1503.28 

30.70 

1511.99 

1511.37 

1507.67 



2ZJ1 

1514.18 

1513.88 


1509.89 


34.76 

1516.46 

1516.94 

1512.51 



36.79 

1512.90 

1516.52 


1511.55 

1509.2 

38.73 

1509.48 

1511.82 

1510.93 



40,61 

1507.30 

1510.22 


1509.89 


42.50 

1508.13 

1509.30 

1509.53 



44.47 


1509.39 


1510,17 

1509.89 

46.40 

1510.29 


1509.62 



48.36 

1511.31 

1510.05 


1509.84 


50.27 

1511.36 


1508,78 



52,20 

1511,03 

1509.79 


1509.36 


54,12 

1510.74 


1509,20 



56,03 

1510.92 

1509.99 


1510.66 

1506.56 

59,88 

1510.31 

1509.19 


1510.64 

1506.78 

63,67 

1509.11 

1507.65 


1508.95 

1505.43 

65.62 

1508.11 


1507.86 



67.54 

1509.14 

1507,55 


1507.92 


71.42 

1501.32 

1509,26 


1508.42 

1505.01 





Speed of Sound versus Depth 
for various frequencies (in kHz) 
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the non-bubble case. After doing this it is possible to compare the 
measured values of c with the values for bubble free sea water. The 
plot of the temperature-caused changes has to be referenced to a value 
of c that holds for the bubble-free medium and the previous calculated 
value of 1506.27 m/sec at 14.3 m was used. The non-bubble theoretical 
value of c is subtracted from the measured value and the difference is 
then called ACg, due to bubbles, which is plotted on the same graph. 

Since the shift in the frequency peaks causes the pattern of the measured 
c versus depth plot to be markedly different for the various frequencies, 
five different frequencies have been picked that are felt to be repre¬ 
sentative for the total frequency region: 36.79, 44.47, 56.03, 63.67, 
71,42 kHz. 

This representation of c versus depth at a constant frequency can 
only be used to examine whether one encounters, at a specific depth, 
resonant bubbles such that for sound frequencies smaller than the bubble 
resonance frequency one experiences a decrease of sound speed and for 
impinging sound of higher frequencies than the resonance frequency an 
increase in c. It should be clear that this plot ignores depth effects 
on an identifiable surface-generated bubble carried to different depths. 

Figure 10 presents this information for a frequency of 36.79 kHz. 
Since there is a strong positive ACg of 5.25 m/sec at 4.3 m this means 
that we encounter strong populations of bubbles with lower resonance 
frequencies than the impinging sound, in other words, there are many 
bubbles that resonate below 36.79 kHz. Figure 11 represents the situ¬ 
ation at 44.47 kHz and shows that we sense bubbles of resonance frequency 
greater than the incident frequency at 4.3 m depth and of resonance 
frequency smaller than the incident frequency at 14,3 m. The picture 
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changes when going to Fig. 12. At a frequency of 56.03 kHz there is 
no significant dispersive behaviour at all . For frequency 63.67 kHz, 
Fig. 13, at all depths we sense bubbles of resonance frequencies greater 
than the frequency of the acoustical signal. Since the difference 
between the measured value and the non-bubble theoretical value is 
much larger at 4.3 than at 14.3 m we are either very close in frequency 
to resonant bubbles at 4.3 m or the population is considerably stronger 
than at 14.3 m. Figure 14, 71.43 kHz, finally shows a smaller negative 
ACg with only minor differences between the various depths. 

Summarizing, we can make the following statements related to the 
different depths: 

4.3 m : There is a large population of bubbles having a resonance 
frequency close to and somewhatsmaller than 36.79 kHz. At 44.47 and at 
56.03 kHz we sense and at 63.67 kHz we detect another population of 
bubbles with an greater than 63.67 kHz. 

7.3 m : Bubbles with a resonance frequency less than 36.79 kHz are also 
detected at this depth. At 44.47 kHz we are far enough away from this 
region and encounter no dispersive characteristic. 

9.3 m : There is a population of bubbles resonating below 36.79 kHz. 

14.3 m : There are bubbles of resonance frequencies below 44.47 kHz. 

At 56.03 kHz this region of resonance frequencies has been passed and 
no dispersive behaviour is shown. 

Conclusions : 

At 4,3, 7.3, and 9.3 m we observed bubbles of resonance frequency below 
36.79 kHz and above 71.42 kHz and at 14.3 m we found a broad region of 
bubbles resonating below 44.47 kHz. The sound speed dispersion is 
strongest at 4,3 m and is reauced in amplitude when going to a greater 











depth. The conclusions derived from the variation in the speed of sound 
with depth are consistent with the observations made for the variation 
in the speed of sound with frequency. 




IX. STATISTICAL AND SPECTRAL ANALYSIS 


The signal received at hydrophone one does not have a constant phase 
difference with respect to that at hydrophone two but is a random signal. 
This is the equivalent of a noisy signal where the noise in this case is 
due to the integrated effects of time-varying speed of sound due to 
varying temperature, salinity and bubble patches along the sound path. 

It is assumed that the random signal exhibits statistical stationarity 
which is partially described in terms of a probability density function 
(P.D.F.) and the statistical moments. In addition to a statistical 
description it is valuable to use the time series to determine the 
temporal correlation function and the power spectral densities of the 
time-varying phase (Ref. '2). 

The variable phase signal that we received at a specific frequency 
and for a period of approximately 90 sec is a continuously changing 
random variable within some finite interval. But since this recorded 
information is digitized, the random variable can assume only a finite 
number of distinct values so that we must treat it as a discrete random 
variable. In practice the recording technique used resulted in magni¬ 
tude increments of .254° spaced in time by .04064 sec. For one frequency 
record we then get n data points. If A is one of the possible data points, 
then n^ is the number of times A occurred and the ratio n^/n is the 
relative frequency of occurrence of the event A for that particular 
record. This approach allowed the use of the sampled data points to 
form a histogram that can be interpreted as a distribution function: 

F(x) = PC<}>:^x) 






from which a density function can be obtained: 


p(x) = ^ P(x.) 5(x - X.) 
i 

Examining these histogramsr it turns out that the main pattern of the 
envelope follows a Gaussian distribution. This can be expected since 
in the experimental measurement many irregular and fluctuating causes 
like weak temperature and salinity inhomogeneities caused the phase 
reading to have small random variations which produce the Gaussian 
P.D.F. about the undisturbed value. It seems that the bell-shaped 
curve experiences stronger changes around 35 and 65 kHz. This different 
behaviour is pointed out by a series of historgrams (Figs. 15-19) for 
the latter frequency range and for 4,3 m depth. At 59.9 kHz the general 
shape is narrow Gaussian but at63.7 kHz the distribution is nearly 
evenly spread over a wider center range that sharply drops at the 
skirts. At 67.6 kHz the envelope changes again and achieves at 69.6 
kHz a wide Gaussian pattern. At 77.3 kHz the pattern is again that of 
a narrow Gaussian noise. Since this change in the density function 
is a particular phenomenon for frequencies around 35 and 65 kHz and is 
not found for frequencies in between, it is postulated that time varying 
resonant bubble populations predominantly around frequencies 35 and 65 
kHz caused the variation with frequency. These frequencies are essen¬ 
tially the same as observed previously in the peaks in the analysis of 
speed as a function of frequency. 

The statistical average, or mean value, T of the random variable <}> 
has been used to calculate the mean speed of sound which is the expected 
value for c. By calculating the variance 

2 ~f ~1 - 2 

= (x - x) = X - X 
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Fig. 17. Histogram of Phase Angle 
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and taking the square root of the variance to obtain the standard 
deviation, a , one has a measure of the spread of the observed values. 

A large standard deviation implies a large spread of likely values for 
any given observation, and vice versa. The standard deviation was 
calculated for all frequency readings from 30.7 to 79.2 kHz and is 
plotted as a function of frequency in Fig. 20. 

Since there was a time lag of two minutes between frequency runs, 
a is for this case not only a function of frequency but also of time, 

A 

which is increasing with frequency. To be able to judge which fraction 
of the variation in a is a time effect and which is the frequency 

/N 

effect a continuous 20 min record of phase modulation at 60 kHz was 
split into ten two-minute records so that the standard deviation could 
be examined as a function of time alone (Fig. 21). Comparing the two 
graphs there is a mean value of 1.28° and a total variation of i .32° 
for the time record; there is a much stronger excursion for the frequency 
plot which shows a mean value of 1.92° at 60 kHz with a deviation of 
+ 1.92° at 40.6 and of + 2.56° at 69.6 kHz. One can rule out temperature 
and salinity patches as the causing agency for the effect in Fig. 20 
since those would not show a frequency dependent spread of the measured 
values. It seems that the presence of bubbles in the medium causes this 
strong change in the standard deviation as a function of frequency. The 
frequencies of the peak values of a are approximately the same as the 
frequencies already identified as centers of resonant bubbles. 

Next the autocorrelation function R(t) will be considered since it 
provides information about (!)(t) and also relates to the frequency-domain 
description of the random signal. RCx) is defined by 

R(t) = < x(t)xCt+T)> /<x(t)x(t+T)> 

T=0 
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Fig. 20. Standard Deviation as a Function of Frequency. 
Time increasing with frequency 



Fig, 21. Standard Deviation as a Function of 
Time Only at Frequency 60 kHz 
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It should be noted that the autocorrelation function is a measure of 
both temporal variation and statistical dependence. The plots of the 
normalized autocorrelation at 30.7 and 46.4 kHz are considered to be 
representative and are given in Figs. 22 and 23. The graph for 30.7 
kHz shows a variation superimposed on the almost ideal damped cosine 
pattern' obtained at 46.4 kHz. This variation is probably due to bubble 
effects. Recall that the frequency 30.7 kHz was observed as a predo¬ 
minant resonance frequency in the speed of sound discussion of section 
VIII. Ihe time it takes for the autocorrelation to drop to 1/e of its 
value is plotted versus frequency on Fig, 24, No specific regularity 
can be found. On the average this "correlation time" is approx. 1.8 sec. 

Examining the autocorrelation function we arrive at the following 
results: The phase modulation remains correlated (1/e) for 1.46 sec, 
the salinity for 3.96 sec and the speed of sound as measured by the 
velocimeter for 2.28 sec (Figs. 25-27). The temperature records have 
not yet been analyzed (thesis of LCDR Duchock to be issued March 1972). 
Since for the velocimeter (operating at 3 MHz) c = c (T, S) it appears 
that the temporal variation in c is due principally to salinity patches 
(and most probably temperature patches). On the other hand, our phase 
measurement is a function of not only S and T but also of the bubble 
distribution. The still smaller correlation time of 1.46 sec is evidence 
of additional decorrelation due to bubbles. 

Turning to the frequency domain, the power spectral density of the 
phase variation gives the distribution of average power of the fluctu¬ 
ations in frequency. To find the PSD the Mean Lag Product Method by 
B1ackman/Tukey was used. The main feature of this method is to calculate 
not only the autocovariance function but in addition an apparent 
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Fig. 22. Autocorrelation at 30,7 kHz 
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Fig. 23. Autocorrelation at 46.4 kHz 















o 

00 



CSJ O 00 

CM CSJ I— I— 


CM 


4 - 


[Das]-3WIi 


Fig. 24. Decorrelation Time as a Function of Frequency 
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Fig. 25. Autocorrelation for Phase Modulation at gO 
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Fig. 26, Autocorrelation for Speed of Sound Fluctuations 
by Velocimeter 
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Fig. 27. Autocorrelation for Salinity Fluctuations 





























autocovariance by using a Parzen-lag-window. By taking the Fourier 
transform of this function one gets a smoothed power spectrum. The 
result of this computation using computer program (2) is a plot of phase 
angle squared per Hz versus frequency (The frequency resolution is .05 
Hz. For better resolution in the amplitude, 10 log-jQ of the ordinate 
is taken such that the dimension is in dB). Several conclusions can 
be derived from the PSD-plot for the 4,3 m depth. For nearly all 
frequencies there is a distinct peak around .48 Hz probably due to the 
orbital effects caused by the observed predominant ocean wave component. 
From 0.5 to 2.5 Hz the values for the PSD follow an exponentially 
decaying curve. This pattern is interrupted by peaks mainly at .8, 

1,3, 1,55, 1.95, and 2.35 Hz that may show up at two or at most three 
consecutive readings (see peak at 2.3 Hz in Figs. 28-31), or occur at 
one frequency and the next time two records later. Since every recording 
took two minutes this then suggests that the agency forcing these changes 
in the frequency spectrum is of a nature that can result in variations 
that show up only over a time period of four minutes or repeat every 
four minutes. 

In general temperature, salinity, and/or bubble patches carried by 
turbulence and crossing the acoustical path could be the driving force 
for these variations. Since the sound signal integrates all these 
effects it is necessary to examine the temperature and salinity records. 
For the 20 min run and the depth of 4.3 m the power spectral density 
of the phase modulation with a frequency resolution of .0045 Hz, of the 
Ramsay velocimeter reading, and of the Bissett-Berman salinometer, all 
follow approximately the same exponentially decaying pattern. Figs. 

31-33 (these graphs courtesy of LCDR Ouchock and LT Seymour). The 
principal difference is that the PSD of the phase shows more power 
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Fig. 28, Power Spectral Density, 46.4 kHz 
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Fig. 29, Power Spectral Density, 48.4 kHz 
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Ftg. 30. Pov/er Spectral Density, 50.3 kHz 
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Ftg. 31. Power Spectral Density for 
Phase Modulation at 60 LHz 
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Fig. 32. Power Spectral Density for Speed of Sound 
Fluctuations by Velocimeter 
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Ftg. 33. Power Spectral Density for 
Salinity Fluctuations 



































concentrated below 0.5 Hz where orbital velocities are strong. All 


three plots of 
the range from 


PSD show good agreement even in the fine structure i 
0.5 - 1.0 Hz. 










X. ERROR ANALYSIS 


The speed of sound as derived before is given by 

c = lx = fx 

n ± (i>/360 

In order to define the accuracy with which c is given we use the 
following approach: 

log c = log f + log x - log Cn ± (i)/360) 

dc = ^ ^ „ d(n ± (i)/360) 

^ f X ^ ^ 

then the maximum error is given by 

Ac = ^ f M + A(n ± (i)/360 

^ f X ^ ^ 

This error will be calculated for the highest (80 kHz) and the lowest 
(30 kHz) frequencies. The coherent decade frequency synthesizer 
synthesizes the frequency of the output signal by proper combination 
of a number of frequencies, each derived on a proportional basis from 
a single internal master frequency crystal. It was essential in this 
application because of the requirement for a stable, sine-wave signal 
source capable of ultraprecise frequency adjustments through a wide 
frequency range. Since the output is always coherent in frequency with 
that of the quartz-crystal oscillator source it has the same percentage 
accuracy as the source which is essentially the smallest digit in the 

_7 

specific frequency or 1.6 x 10 at 60 kHz. This then implies that 

-2 

= 3.3 X 10'^ at 30 kHz and 

30x10"^ 
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at 80 kHz 


10 


-2 


80x10' 


= 1.25 X 10 


-7 


Next the Ax/x term will be examined. The absolute hydrophone distance 
could only be measured to one part in 1.5 x 10 (0.5 mm) which is 

equivalent to 


AX 

X 


5x10 


-4 


= 6.4 X 10 


-4 


,78 


Assuming the Af/f and the At|j/i|) terms to be zero this error in measuring 
the separation between the hydrophones results in an absolute error for 
the speed of sound of 1 m/sec (a level shift of 1 m/sec in the non¬ 
bubbly value of c), but does not affect the dispersion curve, that is 
the pattern of the speed of sound versus frequency. Therefore this 
error is no serious problem and will be ignored in the analysis since 
we are not interested in the absolute accuracy but rather the dispersion 
in the speed of sound. 

To reduce possible errors due to changes in the geometrical config^ 
uration, especially in the separation distance between the hydrophones, 
it was necessary to use as a mounting a wire put under 100 lb of tension. 
The remaining possible vibration due to forces at sea was estimated to 
be 0.1 mm and results in 


AX 
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1.28 X 10“^ 


The accuracy of the phase was the limiting factor. Since basically 
a phase interference approach was used to measure the speed of sound 
any relative phase changes could not be allowed to occur in the two arms 
of the system. It is a question whether one can assume the two parallel 






parts within the system (hydrophone-preampiifier-filter-input channel 
phase meter) to be identical with respect to their influence on the 
signal. Actually, complete identity can never be obtained over a wide 
frequency range. Therefore special attention was given to the phase 
characteristic in selecting the hydrophones. The LC-5 resp. 

LC-10 were selected because their resonance is far above the frequency 
range of interest. Comparing the dispersion curve for 4.3 m with that 
at 14.3 m also points out that there is no phase variation in the 
pattern that shows up consistently at all depths and may therefore be 
due to a different phase behaviour of the two hydrophones. In addition 
the zero control of the phase meter has been used to compensate at a 
single frequency for unequal phase in the two receiving parts when the 
two input signals were known to be exactly in phase by using the 
electrical signal from the output of the frequency synthesizer. 

The accuracy of the analog output of the phase meter is ±0.1°±0.3% 
of the phase angle measured for all frequencies or ±0.16° for the 20° 
phase angle which was used normally in the experiment. This causes a 
fractional error 


a(n ± (|)/360) ^ .16 _ 

n ± <|>/360 ' (360 x 16) + 20 


= .272 X 10~^ at 30 kHz 


.16 


(360 X 41) + 20 


.108 X 10“^ at 80 kHz 


There are of course accuracy limitations in measuring the phase of 
distorted wave forms. Because the phase meter makes its phase measure¬ 
ments by measuring the fraction of a cycle between the negative-going 
zero axis crossing of the two test signals harmonic or scattered signal 









can displace the zero crossing in either a plus or a minus manner 
depending on the phase relationship to the test signal. 

The scattered signal was measured to be dov/n by 37.4 dB re 1 V/y bar 
as compared with the main signal which for the worst case could result 
in an error of .77°, A scattered signal a' adding to the main signal b 
would result in the 1 argesterror in e when the condition depicted 
below occurs 


such that 

for small 0 
then 

and 

A possible source of error is also the fact that at some frequencies 
the response of one hydrophone is significantly stronger than the other 
which means the possibility of harmonic content in the signal at 
multiples or submultiples of these frequencies. The limit error for 
a harmonic of magnitude, e, is 0 = — 57° where a is the magnitude of 

d 

the fundamental (manual for digital input phase meter model 355). 

Since the second harmonic was measured to be -46.34 dB re 1 V/ybar 
as compared with the fundamental the harmonic content could at most 
result in an error of 0.28°. 

Assuming the worst condition for the phase reading or an error of 
.16° at the phase meter, a maximum additional error of ,77° due to the 



0 = - .0135 in radians 

t ,77° 






scattered signal, and also a maximum error of .28° due to harmonic 
content in the signal the total phase error could be 


A(n ± 4>/360) 

n + <t)/360 


.16 + ,77 + .28 

(360 X 16) + 20 


2.08 X 10‘^ at 30 kHz 


1.207 _ 

(360 X 41 ) + 20 


.82 X 10“^ at 80 kHz 


Summarizing it can be stated 

— = 3.3 X 10“^ + 1.28 X 10'^ + 2.08 x 10“^ = 3.36 x 10~^ at 30 kHz 

c 

= 1.25>xl0"^ + 1.28 X 10~^ + .82 x 10“^ = 2.1 x 10“^ at 80 kHz 


Therefore the maximum error in c is mainly determined by the percentage 
error in x and in <}) and is given by 

Ac = (1500) (3.36 X 10^^) = .5 m/sec at 30 kHz and 

= (1500) ( 2.1 X 10“^) = .31 m/sec at 80 kHz 

Comparing the value for Ac at the lov/er and at the upper limit points 
out that in general one is able to take more accurate measurements at 
the higher frequencies. 









XI. SUMMARY AND CONCLUSIONS 


A fixed source-hydrophone system separated by 1.8 m was used to 
study a) the in-situ speed of sound as a function of frequency and 
depth and b) the temporal characteristics of the phase fluctuations 
of a continuous wave acoustical signal over the frequency range 25 - 80 
kHz. The experiment was carried out during Sea State 1 conditions at 
the NUC Oceanographic Research Tower in water of depth 56 ft one mile 
off shore at Mission Bay/San Diego in October 1971. The near surface 
region is of particular interest because bubbles and turbulence can 
have significant acoustical effects in addition to those caused by 
temperature and salinity inhomogeneities. 

Four different depths, 4.3, 7.3, 9.3, and 14.3 m, were examined. 

The near surface values of c ranged from + 6 m/sec to -3 m/sec relative 
to the value for bubble-free sea water, with a maximum error of 0.5 m/sec. 
The speed of sound was found to peak at approximately 35 and 71 kHz at 
4.3 m depth. The lower peak frequency increased as the depth is in¬ 
creased but not at the rate which would be expected for bubbles that 
are wind-initiated at the surface and carried to greater depths without 
losing their identity. A possible explanation is that the bubbles 
originated within the volume or at the bottom and lost gas by diffusion 
as they rose to the surface. Another possibility is that presence of 
more high frequency bubbles at the surface affects the apparent 
resonance frequency of the larger bubbles as determined by sound dis¬ 
persion. The gradual and smooth drop-off in the resonant behaviour with 
increasing depth demonstrates that bubbles of this average size are 
somewhat continuously distributed in the medium down to 14.3 m. The 
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second distinct peak centered at 71 kHz appears to be only a character¬ 
istic of the shallov/est run and does not repeat at the other depths. 

It is postulated that these small bubbles are introduced at the surface 
either by randomly breaking waves or by continental aerosols that drop 
into the water and carry tiny air bubbles with them (Medwin, 1970). 

The surface initiated bubbles disappear quite fast due to gas diffusion 
out of the bubble when transported downwards in the medium. We conclude 
from the foregoing that we encountered bubbles of radii centered around 
124 and 54 microns and resonating at 30 and 69 kHz respectively. The 
broad dispersive characteristic around 30 kHz implies many bubbles 
larger than 60 microns whose principal source was at the surface and 
which were entrained by convective currents carrying them downwards. 

Statistical analysis shows that for frequencies other than at the 
peak bubble populations the phase fluctuations can be approximated by 
a Gaussian probability density function. However around the peak 
variations in sound speed the PDF loses the smooth Gaussian appearance 
and shows a wide, somewhat irregular spread. The temporal variation 
of the standard deviation of the sound phase was small (t .32°). How¬ 
ever, the standard deviation of the phase fluctuation changed 
considerably as a function of frequency increasing to 1.92° at 40.6 
kHz and to 2.56° at 70 kHz, the two sound dispersion centers of this 
experiment. 

Typically, the phase was decorrelated Cl/e) after 1.46 sec. This 
compares with preliminary oceanographic information which shows corre¬ 
lation times of 3.96 sec for the salinity and 2.28 sec for the speed 
of sound as measured by the velocimeter. These values show internal 
consistency. 













♦ 

.. 






The power spectral density of the varying phase shows its strongest 
values at frequencies less than 0.5 Hz; this is presumably due to 
bubbles entrained at orbital frequencies associated with the surface 
wave system which has maximum energy in the same frequency range. 

Due to time considerations it was not possible to further examine 
the phase fluctuations in terms of the functional dependence on 
temperature, salinity, turbulence, and wave height. The analysis of 
the oceanographic parameters will be carried on in the theses by 
Lt. Bordy, LCDR Duchock, and Lt. Seymour (March 1972) and the emphasis 
should then be placed on the interrelation of the parameters as given 
by the cross-correlations or cross spectra for instance. 


{ 
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XII. RECOMMENDATIONS 


It is recormiended that research in this area be continued, especially 
by making an in-situ determination of the bubble concentrations at the 
same time as sound speed measurements are taken in order to determine 
the variation of bubble population with location, v/eather conditions, 
day or night, and local effects. It is also suggested to examine a 
smaller frequency range in more detail rather than to skim through a 
wide frequency range in larger steps and thereby losing valuable 
information. Frequencies below 30 kHz seem to be most inviting because 
of relevance to the Navy and promising because the largest changes in 
the speed of sound occur here, just below the frequency of the most 
prominent bubble population. 








APPENDIX A. PRESENTATION OF DATA 


Run 6, 4.3 m 


freq 

[kHz] 

real t 
[sec] 

mean 

[Dig.] 

standard 

dev. 

[Dig.] 

mean 

[radians] 

n 

fx 

c=-5^ x=.783 

30.72 

180 

-.84 

.06 

-.07467 

15.0987 

1511.99 

32.70 

101 

-.83 

.08 

-.07378 

16.9096 

1514.18 

34.71 

90 

-.69 

.06 

-.06133 

17.9220 

1516.46 

36.82 

80 

+ .445 

.11 

+.03956 

19.0562 

1512.90 

38.70 

30 

+ .65 

.08 

+.05778 

20.0745 

1509.48 

40,60 

38 

+ .83 

.12 

+.07378 

21.0905 

1507.30 

42.42 

44 

+ .08 

.08 

+.00711 

22.0238 

1508.13 

46.40 

193 

+ .44 

.06 

+.03911 

24.0558 

1510.29 

48.40 

107 

+ .665 

.06 

+.05911 

25.0758 

1511 .31 

50.30 

102 

+ ,48 

.06 

+.04267 

26.0593 

1511.36 

52.22 

96 

+ .485 

.03 

+.04311 

27.0598 

1511.03 

54.14 

102 

+ .49 

.06 

+.04356 

28.0602 

1510.74 

56.09 

99 

+ .57 

,06 

+,05067 

29.0673 

1510.92 

58.00 

98 

+ .49 

.06 

+,04356 

30.0602 

1510.77 

59.92 

90 

+ .54 

.06 

+.04800 

31.0647 

1510.31 

61.82 

89 

+ .53 

.06 

+.04711 

32.0638 

1509.65 

63.72 

89 

+ .50 

,07 

+,04444 

33.0611 

1509.11 

65,64 

103 

+ .51 

.11 

+.04533 

34.0620 

1508.90 

67.59 

121 

+ .58 

.11 

+.05156 

35.0682 

1509.14 

69.60 

104 

+ .66 

.14 

+,05867 

36.0753 

1510.64 

71.56 

86 

+ .65 

,10 

+,05778 

37.0745 

1511.32 

73,49 

90 

+ .70 

,11 

+,06222 

38.0789 

1511.14 

75,38 

101 

+ .59 

.09 

+,05244 

39.0691 

1510.72 

77,28 

no 

+. 56 

.07 

+.04978 

40.0665 

1510.25 

79.20 

78 

+ .57 

.07 

+.05067 

41.0673 

1510.05 



J 



Run 7, 

4.3 m 






freq 

[kHz] 

real t 
[sec] 

mean 

[Dig.] 

standard 

dev, 

[Dig.] 

mean 

[radians] 

n 

f X 

c=— x=.783 
n 

30.69 

73.6 

-.94 

.06 

-.0836 

15.8997 

1511.37 

32.69 

112.8 

-.85 

.04 

-.0756 

16.9077 

1513.88 

34.69 

113.6 

-.87 

.05 

-.0773 

17.9060 

1516.94 

36.60 

88.0 

-.97 

,04 

-.0862 

18.8971 

1516.52 

38.76 

98.4 

+ .65 

.09 

+.0578 

20.0745 

1511.82 

40.61 

93.6 

+ .43 

.07 

+.0382 

21.0549 

1510.22 

40.61 

93.6 

+ .44 

,07 

+.0391 

21.0558 

1510.16 

42,53 

104.8 

+ .53 

.05 

+.0471 

22.0638 

1509.30 

44.45 

93.6 

+ .47 

.04 

+.0418 

23.0585 

1509.39 

48.34 

88.8 

+ .55 

.05 

+.0489 

25.0656 

1510.05 

52.19 

92.8 

+ .56 

,06 

+.0498 

27.2667 

1509.79 

56,04 

104.8 

+ .48 

.06 

+.0427 

29.0594 

1509.99 

59.86 

78.4 

+ .45 

.09 

+.0400 

31.0567 

1509.19 

63.66 

99.2 

+ .51 

.09 

+.0453 

33.0620 

1507.65 

67.50 

99.2 

+ .47 

,07 

+.0418 

35.0585 

1507.55 

71,44 

73.6 

+ .52 

.09 

+.0462 

37.0629 

1509.26 

75.24 

92.0 

+ .50 

.08 

+.0444 

39.0611 

1508.22 

79.10 

105,6 

+ .51 

.07 

+.0453 

41.0620 

1508.34 









Run 8, 9.3 m 


freq 

[kHz] 

real t 
[sec] 

mean 

[Dig.] 

standard 

dev 

[Dig.] 

mean 

[radians] 

n 

c=— x=.783 
n 

25.17 

81.6 

+ .51 

.03 

+.0453 

13.062 

1508.81 

28.98 

79.2 

+ .60 

.03 

+.0533 

15.070 

1505.73 

32.91 

84.8 

+ .56 

.05 

+.0498 

17.0665 

1509.89 

36.79 

96,0 

+ .46 

.03 

+.0409 

19.0576 

1511.55 

40.62 

91.2 

+ .54 

,03 

+.0480 

21.0647 

1509.89 

44.49 

93.6 

+ .57 

.06 

+.0507 

23.0674 

1510.17 

48.34 

102.4 

+ .59 

.07 

+,0524 

25.0691 

1509.84 

52,18 

95.2 

+ .59 

.04 

+.0524 

27.0691 

1509.36 

56.07 

103.2 

+ .51 

,10 

+,0453 

29.0620 

1510.66 

59.94 

88.8 

+ .58 

.08 

+.0516 

31.0683 

1510.64 

63.71 

101.6 

+ .48 

.08 

+.0427 

33.0594 

1508.95 

67,56 

108.0 

+ .49 

.07 

+.0436 

35.0603 

1507.92 

71.39 

54.4 

+ .46 

.05 

+.0409 

37.0576 

1508.42 







Run 9, 

14.3 m 



freq 

[kHz] 

real t 
[sec] 

mean 

[Dtg.J 

standard 
dev, 
[Dig.] 

28.95 

42.4 

+ .70 

.05 

36.76 

28.0 

+ .62 

,08 

44.48 

32.0 

+ .56 

.04 

55.93 

36.8 

+ .58 

.06 

59,78 

36.8 

+ .54 

.05 

63.57 

25.6 

+ .53 

,06 

71,28 

78.4 

+ .76 

.09 


mean 

[radians] 

n 

x=.783 

+.0622 

15.0789 

1503.28 

+.0551 

19.0718 

1509.20 

+.0498 

23.0665 

1509.89 

+ .0516 

29.0683 

1506.56 

+.0480 

31.0647 

1506.78 

+.0471 

33.0638 

1505.43 

+.0676 

37.0843 

1505.01 






Run 10, 

7.3 m 



freq 

[kHz] 

real t 
[sec] 

mean 

[Dig.] 

standard 

dev, 

iDig.J 

27,05 

81.6 

.59 

.04 

30.96 

71.2 

.70 

.05 

34.89 

96,8 

.51 

.04 

38,72 

89.6 

.55 

.04 

42,55 

87.2 

.61 

.04 

46.40 

89.6 

.56 

,05 

50.23 

81.6 

.57 

,09 

54,09 

86,4 

.52 

,05 

57.70 

84,8 

'^l ,08 

.07 

61.59 

80,0 

-.79 

.06 

65,59 

78,4 

00 

..13 


mean 

[radians] 

n 

x=.783 

.0524 

14.0691 

1505.44 

.0622 

16.0789 

1507.67 

.0453 

18.0620 

1512.51 

.0489 

20.0656 

1510.93 

,0542 

22.0709 

1509.53 

.0498 

24.0665 

1509.62 

.0507 

26.0674 

1508.78 

.0462 

28.0629 

1509.20 

-.0960 

29.8873 

1511.65 

-.0702 

31.9131 

1511.13 

.0427 

34.0594 

1507.86 
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PROGRAM 1 

PROGRAM CONVERKINPUT»OUTPUT »PUNCH) 

AVGX IS THE AVERAGE DISTANCE DISTANCE THE RECORDING 
DRUM ADVANCES PER ONE HOUR REAL TIME, I.E. THE DATA 
SAMPLING RATE, 

DIMENSION U(6D00»,V(6000),N(80),IBUFFI5000),NK(80) 
IK=0 

READ 97,L 

97 FORMAT I 12) 

READ 98,AVGX 

98 FORMAT(FIO.O) 

READ 96,IFILE 

96 FORMAT!12) 

DO 116 JJ=1,L 
DELT=36.0/AVGX 
PRINT 99,DELT 

99 FORMAT(1H0,2X,25H SAMPLING INTERVAL EQUALS, 1X,E15.7, 
11X,7HSEC0NDS) 

PRINT 700,JJ 

700 F0RMAT(1H0,9HJJ EQUALS, 2X,I2) 

DO 100 1=1,5500 
U(I)=0.0 

100 V(I)=0.0 

101 DO 202 1=1,5000 
202 IBUFF(I)=0 

COUNTX=0.0 

COUNTY=0.0 

NUM=5000 

M=1 

K8=0 

CALL LIOFI5LRBC01,IBUFF,NUM,NPAR,NEOF) 

IF(NEOF) 602,200 

200 K=-7 
KF=0 

201 K=K+8 
KA=K+8 
KC=0 

DO 102 KB=K,KA 
1/ r — 1/ r 4.1 

Nk7kC)=IBUFF(KB) 

102 IF INK(KC).EQ.0) KF=1 
DECODE(80,103,NK) IN(I),1=1,80) 

103 FORMATiaORl) 

DO 104 19=1,80 

IFINI19).EQ.50B) GO TO 106 
IFINI19).EQ.55B) GO TO 107 
IFINI19).E0.34B) GO TO 108 

SYMBOL / 150) REPRESENTS AN INCREMENT TRAVEL IN THE 
MINUS X OR Y OIPxECTION BY THE DIGITIZER, 

SYMBOL 0,I55B), REPRESENTS A ZERO INCREMENT TRAVEL 
IN THE X OR Y DIRECTION BY THE DIGITIZER 
SYMBOL 1,1345), REPRESENTS AND INCREMENT TRAVEL IN 
THE POSITIVE X OR Y DIRECTION BY THE DIGITIZER 
SYMBOL *,I47B), IS A FLAG INSERTED IN IN THE RECORD 
BY THE PERSON DIGITIZING BY USE OF THE IRG 
GO TO 104 

106 RX=-0.01 
K8=K8+1 
GO TO 109 

107 RX=0.0 
K8=K8+1 
GO TO 109 

108 RX=0.01 
K8=K8+1 

109 K3=K8/2 
K3 = 2=!'K3 

IFIK3.EQ.K8) GO TO 111 
COUNTX=COUNTX+RX 
IFICOUNTX.NF.0.0) GO TO 110 
GO TO 104 

110 UIM)=COUNTX 
GO TO 104 




oo 


111 COUNTY=COUNTY+RX 
IF(COUNTX.NE.O.O) GO TO 112 
GO TO 104 

112 V(M)=C0UNTY 
COUNTX=0.0 
M=M+1 

IF(M.GT.6000) GO TO 600 

104 CONTINUE 
IFIKF.EQ.l) GO TO 113 
GO TO 201 

113 MAX=M-7 

105 PRINT 205tM,I9 
205 FORMAT(1H0,2I10) 

TOTAL TIME OF THE RECORD EQUALS NUMBER OF DATA 
POINTS TIMES THE SAMPLING INTERVAL, DELT. 
TIME=(M«DELT)/3600.0 
PRINT llStTIME 

115 FORMAT!IH,IOXtZBH TOTAL TIME OF RECORD EQUALS,2X, 
1E15.7,2X,7H HOURS.) 

PRINT 215,(V(I),I=1,M) 

PRINT 215,(U(I),I=1,MJ 
215 FORMATdH ,10X,14F7.2) 

PUNCH 213,(V(I),I=1,M) 

213 FORMAT!14F5.2) 

PUNCH 2100 

2100 FORMAT ( 20H**********=«'^**=f'^=«'*=!'=<') 

116 CONTINUE 
STOP 

600 PRINT 601 

601 FORMAT! 1 HO,20X ,37H=i'*#=»=* U AND V SPACE INADEQUATE 

602 IK=IK+1 

IF!IK.GT.I FILE J603,200 

603 STOP 
END 




no o oooooooooooooooooooo 


PROGRAM 2 

PROGRAM FOR THE SPECTRAL ANALYSIS OF APERIODIC RECORDS 
DATA MUST BE SAMPLED AT EQUAL INTERVALS OF DT 
CALCULATES EITHER THE AUTO SPECTRA OR THE CROSSSPECTRA 

IF NFLAG1=0 t AUTOSPECTRA.IF NFLAG1 = 1 COSPECTRA 

IF NFLAGP = 1 . CALL PRESS, TOTE 

IF NFLAGC = 1 .CALL URMS 

NTS= NUMBER OF TIME SAMPLES 
MLAG = NUMBER OF TIME LAGS 
DT = TIME INCREMENT 
DHZ = FREQUECY SPACING 
FBHZ = LOWEST FREQUENCY (HZ) 

FEHZ = HIGHEST FREQUENCY OF INTEREST (HZ) 

TM =DT=«=MLAG 
DHZ = 1/(2*A*TM) 

FEHZ = B/(2=!=TM)= B/( 2«MLAG=f=DT) 

FEHZ = A*B*DHZ 

A=«'B = TOTAL NO. FREQUENCY SPACINGS OF ENERGY SPECTRUM 
A AND B ARE INTEGER CONSTANTS 
TRANS— INSTRUMENT AND REMARKS 

REAL*8 ITITLE(12),ITITEL(12),ITITIL(12) 

DIMENSION Fl(6000),F2(6000) 

DIMENSION PHI(600),TAU(600),SPHI(600),APHI(600), 

1PHN(600) 

DIMENSION FREQ(700),CYCL(700),SPEC(700),SSP(700) 
DIMENSION CSPEC(700),QSPEC(700),SS2(700) 

DIMENSION RESPF(700),PER(700),SPE2(700).PHASE(700), 
1COHER(700) 

DIMENSION DATE(2),H0UR(2),TRANS(10) 

REAL*4 LABEL/4H CO /,LA8LI/4HQUAD/,MONT,LABIL/4H PHI/, 
1LABLE/4HCR0S/ 

READ(5,200,END=888) (ITITLE(I),1=1,12) 
READ(5,200,END=888) (ITITEL(I),1=1,12) 

200 FORMAT(6A8) 

READ(5,802) DATE,HOUR,TRANS 

802 F0RMAT(2A4,2A4,10A4) 

READ(5,96) NFLAGl,NFLAGP,NFLAGC 
READ!5,99) NTS,MLAG,DT,DHZ,FBHZ,FEHZ 
READ(5,801) CALX1,CALX2,H,X2,AZZ,S 
WRITE(6,811) 

811 FORMAT(lHl) 

WRITE(6,803) DATE,HOUR,TRANS 

WRITE!6,98) NTS,MLAG,DT,DHZ,FBHZ,FEHZ,CALXl,CALX2,H,X2 
96 FORMAT (316) 

98 F0RMAT(7H NTS = ,I6//,8H MLAG = ,I6//,6H DT = ,F9.5//, 
16H DHZ= ,F10 

2.8//,' FBHZ = '.FIO.S//,' FEHZ = *,F12.8//,* CALXl = • 
3,F10.8//, 

4* CALX2 = '.FIO.B//,* H = *,F10.3//,' X2 = '.FIO.B//) 

99 FORMAT (2I5,5F10.8) 

801 F0RMAT(6F10.5) 

803 F0RMAT(40H SOUND PHASE STUDY OF SAN DIEGO BAY , 

15X,7H DATE ,2A4,10H HOUR ,2A4//,5X,10A4///) 

BAND WIDTH FREQUENCIES OF TOTAL ENERGY FLUX- CM1N,CMAX 
CMIN = 0.0 
CMAX = 0.2 

COMPUTING POWER SPECTRUM 
PI=3.14159265 
FB = FBHZ-2.0*PI 
FE = FEHZ-2.0*PI 
DF = DHZ *2.0*PI 
FMIN = 2.0*PI*CMIN 
FMAX = 2.0-PI*CMAX 
900 FORMAT!14F5.2) 

READ(5,900)(FI(I),I=1,NTS) 

XMAX = FKl) 

DO 23 1=2,NTS 

IF(XMAX.GE.f^l( I ) ) GO TO 23 
XMAX = FI(I) 

23 CONTINUE 

DO 24 1=1,NTS 








Fl( I) = (F1(I)-XMAX) 

24 CONTINUE 

WRITE(6,901)(Fl(I),1=1,NTS) 

901 FORMAT!14F7.2) 

CALL TREND(F1 ,NTS,DT,CALX1) 

21 DO 10 M=1,MLAG 
SUM=0.0 
NMAX=NTS-M+1 
DO 8 I=1,NMAX 
NN=M+I-1 

8 SUM=SUM+F1( I)=«=F1(NN) 

XNMAX=NMAX 

XX=M-1 

TAU(M) = XX=!=DT 
PHI(M)=SUM/XNMAX 
PHN<M) = PHI(M)/PHI(1) 

10 CONTINUE 

CALL PARZ(MLAG,PHI) 

798 FORMAT!• PHI!0) = •,F10.4//) 

NFREQ=!FE-FB)/DF+0.1 
DO 14 N=1,NFREQ 
XN=N 

FREQ!N) = !XN-1.0)*DF+FB 

14 CYCL!N) = FREQ!N)/12.0=!'PI) 

PER!1) = 0.0 

DO 15 N = 2,NFREQ 

15 PER!N) = 1.0/CYCL!N) 

MLAGM1=MLAG-1 

XMLAG=MLAG 

46 DO 50 N=1,NFREQ 

SUM=0.5*!PHI!1 )+PHI!MLAG)=«=COS!FREQ!N)*TAU!MLAG) ) ) 
C1=C0S!FREQ!N) ^i'DT) 

S1=SIN!FREQ(N)=!=0T) 

CC=1.0 

SC=0.0 

DO 49 M=2,MLAGM1 

CT=CC*Cl-SC>i'Sl 

ST=SC*C1+CC*S1 

CC=CT 

SC=ST 

49 SUM=SUM+PHI !M) ^‘CC 

50 SPEC!N) = SUM*2. 0/XMLAG 
IF!NFLAGP) 741,741,742 

742 CONTINUE 

CALL PRESS!FREQ,SPEC,NFREQ,H,X2,DF,FB,FMAX) 

741 CONTINUE 

M0=3.14159265/!DF#DT*XMLAG)+0.1 
DC = DF/!2.0*PI) 

SUM = 0.0 
DO 650 N=1,NFREQ 

650 SUM = SUM + SPEC!N) 

WRITE!6,651) SUM 

651 FORMAT !25H VARIANCE OF SPECTRUM//,3X, 

116H VARIANCE = ,F10.5,10H M2 //) 

WRITE!6,106)!TAU!M),PHN!M),M=1,MLAG) 

INUM=FEHZ/DH2 

DO 9988 1=1,I MUM 

SPEC!I)=10*AL0G10!SPEC!I)) 

9988 CONTINUE 

C CALL DRAW TO PLOT SPECTRUM OF SAN DIEGO 

C AS A FUNCTION OF FREQUENCY 

CALL DRAW!200,TAU, PHN,0,0,LABEL,ITITEL,8.0,0.0,2,0,2 
1,0,8,8,1,LAST) 

CALL DRAW!INUM,CYCL,SPEC,0,0,LABEL,ITITLE,0.0,0.0,0,0, 
10,0,8,8,1,LAST) 

105 FORMAT !10X,24H ENERGY SPECTRAL DENSITY///,5X, 

IlOH FREQUENCY,5X,13H 

2 RAW SPECTRUM,3X,18H SMOOTHED SPECTRUM//) 

106 FORMAT !5X , F10.5,6X ,F10.5) 

WRITE!6,105) 

WRITE!6,106)!CYCL!M),SPEC!M) ,M=1,NFREQ) 

888 RETURN 















END 


SUBROUTINE PRESS(FREQ,SPEC,NFREQ,H,X2,DF,FB,FMAX) 

C SUBROUTINE TO CONVERT PRESSURE SPECTRUM TO ENERGY SPECTRUM 
DIMENSION FREQ(NFREQ),SPEC(NFREQ) 

NMAX = (FMAX-FB)/DF+1.0 
DO 11 1=1,NMAX 

C CALCULATE LINEAR WAVE LENGTH BY NEWTONS METHOD 
IF(FREQ(I)-0.00001) 8,8,7 

7 XKHO = FREQ( I )*FREQ( I )=^'H/9.80 
IF(XKH0-6.3) 5,1,1 

1 XKH = XKHO 
GO TO 9 

5 XKH = SQRT(XKHO) 

3 SH = SINH(XKH) 

CH = COSH(XKH) 

EPS = XKHO-XKH«SH/CH 

SLOPE = -XKH/CH*-2-SH/CH 

DXKH = -EPS/SLOPE 

IF(ABS(DXKH/XKH)-0.0001) 9,9,4 

4 XKH = XKH+DXKH 
GO TO 3 

8 RESPF = 1.00 
GO TO 11 

9 XK = XKH/H 

RESPF = C0SH(XK«H)/(C0SH(XK*(H-X2))) 

IF(RESPF-IO.O) 11,11,12 

11 SPEC(I) = RESPF=!'RESPF=<=SPEC{I ) 

12 RETURN 
END 


SUBROUTINE TREND(FX,NTS,DT,CALXX) 

DIMENSION FX(NTS) 

C CALIBRATIONG RECORD 
DO 104 1=1,NTS 
104 FX(I) = FX(I)*CALXX 
C COMPUTING THE LINEAR TREND 
FNTS = NTS 
SUMF = 0.0 
DO 101 1=1,NTS 

101 SUMF = SUMF + FX( I ) 

SUMFl = 0.0 

DO 102 1=1,NTS 
XI = I 

102 SUMFl = SUMFl + XI=t=FX(I) 

XNMl = NTS-1 

XNPl = NTS+1 

XM = (1.0/DT)^(12.0*SUMF1/(FNTS^XNM1=<'XNP1)-6.0*SUMF/ 
1(XNMl^FNTS)) 

B = SUMF/FNTS-XM*XNPl=!=DT/2.0 
FMEAN = SUMF/FNTS 
WRITE(6,9) FMEAN,XM,B 

9 FORMAT (3X,8H MEAN = ,F10.5,3X,9H SLOPE = ,F10.5,3X, 
113H INTERCEPT= ,F10.5//) 

DO 103 1=1,NTS 
XI = I 

103 FX(I) = FX(I) - (B + XM=f'XI*DT) 

RETURN 

END 


SUBROUTINE SMO(MD,XI,X2,NFREQ) 
DIMENSION Xl(MD),X2(MD) 

DO 1 N=1,M0 
NA=N+MD 
NN=NFREQ-N+1 
NB=NN-MD 

X2(N) = 0.25«(X1(1)+X1(NA))+0,5*Xl(N) 
1 X2(NN)=0.5*(X1(NN)+X1(NB)) 

3 MB=MD+1 


I 








ME=NN-1 

5 DO 2 N=MB,ME 
NA=N+MO 
N8=N-MO 

2 X2(N)=0.25*(X1(NA)+X1(NB))+0.5*Xl(N) 
RETURN 
END 


SUBROUTINE PARZ(MLAG,PHI) 

C PARZ SUBROUTINE PARZEN FILTERS AUTO-CORRELLATION FUNCT 

DIMENSION PHKMLAG) 

XMLAG = MLAG 

MLAGH = XMLAG/2.0-0.1 

MLAGHl = MLAGH + 1 

DO 31 M=l,MLAGH 

MM = M-1 

R = MM 

RM = R/XMLAG 

UM = 1.0-6.0*RM*RM*(1.0-RM) 

PHKM ) = PHKM )«UM 

31 CONTINUE 

DO 32 M = MLAGHl,MLAG 
MM = M-1 
R = MM 

RM = R/XMLAG 
RMl = (1.0-RM) 

UM = 2.0=!'RM1*RM1*RM1 
PHKM ) = PHKM )*UM 

32 CONTINUE 
RETURN 
END 







C PROGRAM 3 CALCULATES DISPERSION IN C 

// EXEC FORTCLGP,REGION.G0=100K 
//FORT.SYSIN DD - 

REAL=!'8 TITLE(12)/• BOX 209 Sll*' */ 

REAL LABEL/' •/ 

DIMENSION A(10),B(10),MR(10),W(500),CC(500) 

DATA PI/3.1415926/ 

C SET MAX DIMENSION TO WORK WITH. 

N0UT=500 

C READ PZ,RHZ,CZ. 

REA0(5,100) PZ,RHZ,CZ 
100 F0RMATI3F10.5) 

PZ=PZ=i‘lQOOOO.O 

C READ RADIUS (MICRONS) AND FREQ. 

READ(5,200) ISR,I DR, lER,ISW,IDW,IEW 
200 FORMAT(6I10) 

C READ NO. CF SEGMENTS AND THEIR VALUEX...MAX OF 10 SEGMENTS 
READ(5,300) NS 
300 FORMAT(12) 

REA0(5t400)(A( I),B(I),MR(I),1 = 1,NS) 

400 F0RMAT(2F10.5»110) 

C PRINT INPUT. 

WRITE(6,500) 

500 F0RMAT(•lINPUT:•) 

WRITE(6,510) PZtRHZtCZ 

510 FORMATC PO = * , El 6.8, T28, • RHO = • , E16.8 ,T54,'CO = ', 
lE16.a) 

WRITE(6,520) ISR»IDR,IER 

520 FORMAT!• STARTING RADIUS = ', 16,T28,'INCREMENT = ',13, 
1T54,'ENDING RADIUS = ',16) 

WRITE(6,530) ISW,IDW,IEK 

530 FORMATC STARTING W = ',16 ,T28,*INCREMENT = ',I3,T54, 

1'ENDING W = ',16) 

WRITE!6,540) 

540 FORMATC SEGMENT ' j7X ,'A', 15X, »B» ,9X, ' MAX. RADIUS') 
WRITE!6,550) !I,A!I),B!I),MR!I),1=1,NS) 

5 50 FORMAT!5X,12,3X,El6.8,2X,E16.8,4Xj16) 

C CHECK FOR INPUT ERROR. MR SHOULD INCREASE. 

NJ=NS-1 

IF (NJ .EQ. 0) GO TO 700 
DO 600 J=1,NJ 

IF !MR!J) .LT. MR!J+1)) GO TO 600 
WRITE!6,560) 

560 FORMATC SEGMENT RADII ARE OUT OF SEQUENCE.') 

STOP 

600 CONTINUE 

C NOW CALCULATE SOME CONSTANTS TO SAVE TIME. 

700 IR=(lER-ISR)/IDR+1 
IW=( IEVJ-ISW)/IDW + 1 
WC = 1 .975'-SQRT( PZ/P.HZ)*1000000.0 
PIC=0.027/( (P 1^2000.0) =5^*0.322) 

C CLEAR OUTPUT BUFFER 
DO 800 I=1,N0UT 
W! I )=0. 

800 CC!I)=0. 

C INITIALIZE STARTING "VALUES . 

C WV/=WORKING OMEGA. 

C WZI=CALCULA7ED OMEGA ZERO SUB I. 

C RW=WORKING RADIUS. 

ITEMPW=ISW 
DO 2000 K=1,IW 
WW=ITEMPW 
W!K)=ITEMPW 
IF (K.NE.l) GO TO 915 
WRITE(6,900) ITEMPW 
900 FORMAT(•IW = •,16) 

WRITE(6,910) 

910 FORMAT!' INDEX•,6X,'RADI US',1IX,'WZ',13X,'U(R)D(R)') 

915 SUM=0.0 

ITEMPR=ISR 

C THIS IS THE SUMMATION LOOP. 

DO 1000 L=1,IR 




WZI = WC/ ITEMPR 

DI = PIC=i'(WZI*-0.322 ) 

WRAT=WZI/( WW=i=2 .0=i=P I’i'lOOO.0 ) 

WRAT2=WRAT=^WRAT 
WRATM=WRAT2-1.0 

BATH=WRATM/( (WRATM*WRATM) + (OI’i=DI=<'WRAT2) ) 

C SELECT U(R)D(R). 

DO 920 J=1»NS 

IF (ITEMPR .GT.MR(J)) GO TO 920 
URDR=1.0E-09«(A(J)*ITEMPR+B(J)) 

BUBeLE=URDR/ ( ITEMPR’>=ITEMPR=»:1 .OE-12) 

GO TO 927 
920 CONTINUE 

C IF FALLS THROUGH HEREt WE HAVE GONE TOO FAR. 

925 WRITE(6»926) 

926 FORMAT!* ERROR ON SEGMENT RADIUS.*) 

STOP 

C WRITE RESULT FOR EACH I. 

927 IF (K.NE.l) GO TO 940 
WRITE(6»930) L,ITEMPR,WZI,URDR 

930 FORMAT!2X,I6t6XfI6,6X,E16.8,2X»E16.8I 
940 ITEMPR=ITeMPR+IDR 

SUM=SUM+(BUBBLE*BATH) 

1000 CONTINUE 

C NOW CALCULATE SOUND VELOCITY. 

CNEW=1.0/( (1.0/CZ) + ( (3.0*CZ=<'SUM)/(8.0-*PI*PI*WW^WW* 
11.0E06) ) ) 

CC(K)=CNEW-CZ 

WRITE(6,1010) ITEHPW,SUM,CNEW,CC(K) 

1010 FORMAT!* W = •,I6,4X,*SUM = •,E16.8,2X,•C = *,E16.8, 
12X,'OIFF = •,E16.8 ) 

C PREPARE FOR NEXT ROUND OF W. 

ITEHPW=ITEMPW+IDW 
2000 CONTINUE 

C ALL DONE..NOW PLOT W AND C. 

WRITE!6,2100 ) 

2100 FORMAT!*1 RESULT: W VS. (C-CO)') 

WRITE(6,2200)!W!I),CC(I),I=1,IW) 

2200 FORMAT!10X,E20.10,lOX,E20.10) 

WRITE(6,2300) 

2300 FORMAT!'l') 

CALL PLOTP (W,CC,IW,0) 

STOP 

END 
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